Hydrodynamic slip boundary condition at chemically patterned 
surfaces: A continuum deduction from molecular dynamics 



Tiezheng Qian and Xiao-Ping Wang 
Department of Mathematics, Hong Kong University of Science and Technology, 
Clear Water Bay, Kowloon, Hong Kong, China 

Ping Sheng 

Department of Physics and Institute of Nano Science and Technology, 
Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China 

Abstract 

We investigate the slip boundary condition for single-phase flow past a chem- 
ically patterned surface. Molecular dynamics (MD) simulations show that 
modulation of fluid-solid interaction along a chemically patterned surface in- 
duces a lateral structure in the fluid molecular organization near the surface. 
Consequently, various forces and stresses in the fluid vary along the patterned 
surface. Given the presence of these lateral variations, a general scheme is de- 
veloped to extract hydrodynamic information from MD data. With the help 
of this scheme, the validity of the Navier slip boundary condition is verified 
for the chemically patterned surface, where a local slip length can be defined. 
Based on the MD results, a continuum hydrodynamic model is formulated us- 
ing the Navier-Stokes equation and the Navier boundary condition, with a slip 
length varying along the patterned surface. Steady-state velocity fields from 
continuum calculations are in quantitative agreement with those from MD 
simulations. It is shown that, when the pattern period is sufficiently small, 
the solid surface appears to be homogeneous, with an effective slip length that 
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can be controlled by surface patterning. Such a tunable slip length may have 

important applications in nanofluidics. 
PACS numbers: 83.50.Lh, 83.10.Rs, 83.50.Rp 
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I. INTRODUCTION 



Recently developed techniques of microfluidics [1] have generated great interest in fur- 
ther miniaturization towards nanofluidics. Precise control of fluid flow at nanoscale not only 
challenges existing techniques but also requires more realistic simulations for understand- 
ing nanoscale flow phenomena [2]. Theoretical and experimental studies of confined fluids 
have shown that boundary conditions at various interfaces can strongly influence the flow 
behavior as the system size approaches nanoscale. One particular boundary condition that 
has received substantial efforts is the slip boundary condition at the fluid-solid interface [3] . 
Molecular dynamics (MD) studies of fluid slipping have focused on single-phase flow past 
homogeneous surfaces [4-8] and miscible [9,10] and immiscible [11-13] two-phase flows past 
homogeneous surfaces. On the other hand, chemically patterned surfaces have generated 
interest in novel control of flow in micro- and nanofluidics [14-16]. In particular, slip flow 
past chemically patterned surfaces has been investigated using both MD and continuum 
simulations [17]. Wetting and dewetting on chemically patterned surfaces have also been 
extensively studied [18,19]. 

In this paper we study the slip boundary condition for single-phase flow past a chemically 
patterned surface. MD simulations show that modulation of fluid-solid interaction along the 
patterned surface induces a lateral structure in the fluid molecular organization near the 
surface. As a consequence, various forces and stresses in the fluid vary along the patterned 
surface. Given the presence of these lateral variations, a general scheme is developed to 
extract hydrodynamic information from MD data. With the help of this scheme, the validity 
of the Navier slip boundary condition is verified for the chemically patterned surface, where 
a local slip length can be defined. Continuum calculations are carried out using the Navier- 
Stokes equation and the Navier boundary condition (NBC), with a slip length varying along 
the chemically patterned surface. Changing the ratio of the pattern period to the slip length, 
two distinct slip regimes are found. Particularly, when the pattern period is sufficiently small 
compared to the slip length, the solid surface appears to be approximately homogeneous. 
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with an effective slip length that can be controlled by surface patterning. Such a tunable 
slip length may have important applications in nanofiuidics. 

The paper is organized as follows. In Sec. II, we describe the MD simulation techniques 
for a single-phase fluid confined between two solid planar walls, one with a chemically homo- 
geneous surface and the other a chemically patterned surface. The sohd surfaces are parallel 
to the xy plane and surface patterning is along the x direction (see Fig. 1). Interaction 
between fluid and solid molecules is modeled by a modified Lennard-Jones (LJ) potential 

where e/^ and u/^ are the energy and range parameters, and 5/^ is a dimensionless parameter 
adjusting the wetting property of the fluid. This LJ potential has constant parameters efs, 
afs, and 5fs at the homogeneous surface. For the patterned surface, we choose 5fs to be an 
interaction parameter that varies with the equilibrium x coordinate of the solid molecule 
(see Fig. 1). Denoted by Sfs{x), the variation of 6fs in the x direction models the patterned 
surface in MD simulations. 

In Sec. Ill, we describe how to extract hydrodynamic information from MD data. It 
is shown that in the absence of hydrodynamic motion, modulation of fluid-solid interaction 
along the patterned surface induces a nonuniform static state in which various forces and 
stresses in the fluid vary along the solid surface, and these lateral variations decay quickly 
away from the fluid-solid interface. The forces and stresses measured in the static equilibrium 
state are taken as reference quantities. Together they form a molecular background from 
which hydrodynamic components are generated when the fluid is driven into hydrodynamic 
motion. Couette flow is generated by shearing the fluid. The forces and stresses measured 
in the dynamic steady state are taken as dynamic quantities. A hydrodynamic quantity is 
obtained by subtracting the reference quantity from the corresponding dynamic quantity: 

Q — Qd — Qs: 

where Qs is a static (reference) quantity, Qd the corresponding dynamic quantity, and the 
over tilde denotes the hydrodynamic quantity. 



In Sec. IV, we formulate a continuum hydrodynamic model for slip flow past chemically 
patterned surfaces. Various MD data are collected and analyzed under the subtraction 
scheme described in Sec. III. In particular, the NBC is verified to be a local constitutive 
relation for the hydrodynamic fluid-solid coupling at the patterned surface. According to 
the NBC, the shp velocity is locally proportional to the hydrodynamic tangential stress at 
the wall, and the local proportionality constant defines a local slip length. This may be 
expressed by 

(3{x)vfP{x) = -dnx = -TjdnVxix), 

where v^^^^ is the slip velocity, — 5"„^ = —r]dnVx{x) the hydrodynamic tangential viscous 
stress (?7 is the viscosity and n denotes the outward surface normal) , and 13 [x) the local slip 
coefficient. It is shown that as a local property at the patterned surface, j3{x) is determined 
by the local interaction parameter 5/s(x), as expressed by (5{x) — (3{5fs{x)), where /3 as 
a function of 5fs, denoted by P{Sfs), relates a microscopic interaction parameter and a 
coefficient in continuum boundary condition. The local slip length ls{x), defined by ls{x) = 
T)/ (3{Sfs{x)), varies with the local fluid-solid interaction at the patterned surface. 

The continuum hydrodynamic model formulated in Sec. IV consists of the Navier-Stokes 
equation with a constant viscosity and the NBC with a laterally varying slip length. All 
the parameters involved in the continuum model calculations are directly determined from 
MD simulations, thus no adjustable parameter appears. In Sec. V, we present steady-state 
velocity fields from continuum calculations, in quantitative agreement with those from MD 
simulations. This further affirms the vahdity of the local NBC. Continuum calculations also 
show that there are two distinct regimes for fluid slipping. As the surface pattern period is 
much larger than the slip length, well separated slip domains form along the solid surface. 
When the pattern period is sufficiently small, by contrast, the solid surface appears to be 
homogeneous and the shp velocity can be predicted using an effective slip length. The latter 
regime provides a way to systematically tune the slip length. We conclude this paper in Sec. 
VI. 
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II. MOLECULAR DYNAMICS SIMULATIONS 



A. Geometry and Interactions 

MD simulations have been carried out for Couette flow between two solid walls, one with 
a chemically homogeneous surface and the other a chemically patterned surface (see Fig. 
2). The two planar walls are parallel to the xy plane, with the two fluid-solid boundaries 
defined at 2; = 0, H. The Couette flow is generated by moving the upper wall at a constant 
speed V in the x direction. Periodic boundary conditions are imposed in the x and y 
directions. Interaction between fluid molecules separated by a distance r is modeled by a 
Lennard- Jones (LJ) potential Uff{r) — 4e (cr/r)^^ — (cr/r)^ , where e is the energy scale 
and a the range scale. The average number density for the fluid is set at p = 0.81o"~^. The 
temperature is maintained at T — lAe/ks- Each wall is constructed by two [001] planes 
of an fee lattice, with each wall molecule attached to a lattice site by a harmonic spring. 
The mass of the wall molecule equals that of the fluid molecule m. The number density 
of the wall equals — 1.86cr~^. The fluid-solid interaction is modeled by a modifled LJ 
potential Ufs(r) — 4e/s (afg/ry^ — 5fs (crfg/r)^ , with the energy and range parameters 
given by ejg = 1.16e and afg = l.OAa, and a dimensionless parameter 6fs for adjusting the 
wetting property of the fluid. The lower solid surface is patterned, and the patterning effect 
is modeled by the oscillating parameter [20] 

r Sttx 1 

5fs{x) ^ 5o + 5icos -p- , (1) 

where x is the x-coordinate of the wall molecule, Sq is a constant. Si is the oscillation ampli- 
tude for 5fs, and P is the oscillation period, set to be larger than lOcr. (This requirement will 
be explained in Sec. IV A.) The upper solid surface is homogeneous, at which the fluid-solid 
interaction has a constant Sfg, equaling Sq. We use Sq — 1 and Si = 0.3. The fluid-solid 
interaction potential Ufs is cut off at Tc = 2.2(7. (Thicker walls with larger cut-off distance 
have also been used, but the MD results are quite stable. Sec Appendix A.) In most of our 



simulations, the shearing speed V is on the order of O.l^e/m, the wall separation along z is 
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H — 13.6(7, the sample dimension along y is 6.8(7, and the sample dimension along x is set 
to be L = 2P (from -P to P). The MD results presented in Figs. 3 to 9 and 11 to 14 are 
obtained for ^ = and 0.25^fe/m, H = 13.6a, and L = 2P = 54.5a. 



B. Measurements 



To have spatial resolution along the x and z direcrtions, the measurement region is evenly 
divided into bins, each Ax = 0.85a by Az = 0.85a in size. That is, the fluid space is divided 
along the z direction into H/ Az horizontal layers, and each layer is then divided along the 
X direction into L/ Ax bins. We measure five quantities: G^, the tangential wall force per 
unit area exerted by the wall on the fluid molecules in a horizontal layer; axx and a^xi the 
XX and zx components of the fluid stress tensor; and Vz, the x and z components of 
the fluid velocity. is obtained from the time average of the total tangential wall force 
experienced by the fluid molecules in a bin, divided by the bin area in the xy plane (for more 
details see Appendix B); a^x^zx) is obtained from the time average of the kinetic momentum 
transfer plus the fluid-fluid interaction forces across a constant-x(2;) bin surface (for more 
details see Appendix C); Vx and are measured from the time average of fluid molecules' 
velocities in each bin. As reference quantities, we also measure (7^^., and a^^x ™ 
static iy — 0) state (in which Vx and both vanish). Here the superscript "0" denotes 
the static quantities. Static equilibrium-state {V = 0) and dynamic steady-state {V ^ 0) 
quantities are obtained from time average over 5 x lO^r or longer where r is the atomic time 



scale yma'^/e. Throughout the remainder of this paper, all physical quantities are given in 
terms of the LJ reduced units (defined in terms of e, a, and m). 
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III. STATIC AND DYNAMIC STATES 



A. Static Equilibrium State 

1. Nonuniform static state 

In the present Couette flow geometry, a static equilibrium state is usuaUy characterized 
by zero velocity and constant pressure in a continuum hydrodynamic description. Micro- 
scopically, however, various forces and stresses in the fluid actually vary along the solid 
surface. These lateral variations are induced by the short-range fluid-solid interaction, and 
they decay quickly away from the fluid-solid interface. The periodicity of the surface struc- 
ture necessarily shows up in these near-surface lateral variations. First there is the short 
period of crystalline structure. (For G^", this corresponds to the flrst component discussed 
in Appendix B.) At the lower fluid-solid interface, in addition, there is the large period of 
surface patterning: the fluid-solid interaction itself is periodically modulated along the solid 
wall, with a period much larger than that of the crystalline structure. As a consequence, the 
lateral force/stress variations will show this period of the fluid-solid interaction. (For 
this corresponds to the second component discussed in Appendix B.) In Fig. 3, we show 
the proflles of G^^ for the two boundary layers at the lower and upper fluid-sohd interfaces. 
While G^^ at the upper fluid-solid interface shows only a fast oscillation imposed by the 
crystalline structure, at the lower interface C^^ shows a slow variation as well, induced by 
the periodicity of the fluid-sohd interaction. In Fig. 4, we show the proflle of z-integrated 
for the boundary layer at the lower fluid-solid interface. A slow variation with period 
P is clearly seen. In Fig. 5, we show the proflle of (7°^ z — /^z near the lower fluid-solid 
interface. Once again, a slow variation with period P is clearly seen. 
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2. Static force balance 



In the static equilibrium state = 0), local force balance necessarily requires the stress 
tangential to the fluid-solid interface to be the same on the two sides. Then the wall force 
and fluid stress must vary in such a way that the total force density vanishes. An integrated 
form for the static tangential force balance in the boundary layer at the lower fluid-solid 
interface is given by 



Here l\z/2 represents the mid- level z — /S.z/2 in the boundary layer between z — Q and 
z — Az, G'^^{x, Az/2) is the tangential wall force per unit area in this layer, given by 



in terms of the tangential wall force density g"^^, and Gl^{x, Az/2) is the tangential fluid 
force per unit area in the same layer, given by 



where and a^^ are the xx and zx components of fluid stress in the static state, both 
measured as reference quantities. In Eq. (3) the fact that (7°^(x, 0) — has been used. 
(More strictly, cr^xi^, 0~) = because there is no fluid below z — 0, hence no momentum 
transport across z = 0.) Similar expressions for the static tangential force balance in other 
horizontal layers can be obtained as well. In Fig. 6 we show the profiles of and —Gl'^ 
for the boundary layer at the lower fluid-solid interface, and Eq. (2) is verifled. 

Away from the fluid-sohd interface and deep in the fluid, g'^^ and a^^ vanish while a^^ 
becomes a constant. This leads to another integrated form for the static tangential force 
balance: 



Gf(x,Az/2) + Gi%x,Az/2)^0, 



(2) 




Gi'ix, Az/2) = /o^^ dz[d,al{x, z) + z)] 
= /(f ^ dza^^ix, z) + (7°^(x, Az), 



(3) 



dz[gf{x, z) + d,al,{x, z) + d,al,{x, z)] 

!^'dzgf{x,z)+d,!^'dzal,{x,z) 



= 0, 
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where Zd represents a level z — Zd far away from the interface and deep in the fluid. Therefore, 
for the whole nonuniform fluid layer in equilibrium at the solid surface, there is no net 
tangential force (per unit area), although the fluid-solid interaction is modulated along the 
tangential (x) direction. This is simply due to the cancellation of the wall force by the fluid 
force, because fluid molecules are adaptable to the modulated fluid-solid interaction. 

B. Dynamic Steady State 

1. Hydrodynamic quantities 

In the static state, various forces and stresses show spatial variations in the molecular- 
scale vicinity of the wall. The nonuniform static state forms a microscopic background 
upon which hydrodynamic variations are generated when the fluid is driven (sheared) into 
hydrodynamic motion {V ^ 0). That's why the static quantities are taken as reference 
quantities. 

As observed in MD simulations for 7^ 0, the variation of a particular static force/stress 
is still preserved in the variation of the corresponding dynamic force/stress. To obtain a 
hydrodynamic quantity, we need to subtract the static part from the corresponding dynamic 
quantity. This subtraction scheme is formally expressed as 

Q = Qd- Qs, (4) 

where Qs is a static (reference) quantity, Qd the corresponding dynamic quantity, and 
the over tilde denotes the hydrodynamic quantity. Physically, the characteristic variation 
magnitude of Q should be much smaller than that oiQo and Qs- In Figs. 7 and 8, we show 
the profiles of G'^{x, Az/2), J^"" dza^xix, z), and azxixjAz) for the boundary layer at the 
lower fluid-solid interface. They are obtained according to Eq. (4): 

G:ix, Az/2) = G:ix, Az/2) - Gf{x, Az/2), 

pAz rAz 

/ dzaxx{x,z) = dz 
Jo Jo 
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crxx{x,z) - a^^{x,z) 



and 

It is readily observed that the variation of each hydrodynamic quantity shows a magnitude 
much smaller than that seen from the corresponding static quantity. This fact indicates that 
in deducing a continuum relation (say, a slip boundary condition or a constitutive equation 
for stress), if a hydrodynamic quantity (Q) is to be used, then it can by no means be replaced 
by a dynamic quantity (Qd), without subtracting the static part. This will be shown in 
detail in Sec. IV A for the Navier slip model and Sec. IV B for the Newtonian stress. 

2. Hydrodynamic force balance 

In the dynamic steady state {V 0), local force balance still requires the stress tangential 
to the fluid-solid interface to be the same on the two sides. An integrated form for the 
dynamic tangential force balance in the boundary layer at the lower fluid-solid interface is 

given by 

G:{x,Az/2) + Gi{x,Az/2) = 0. (5) 

Similar to G^^{x, Az/2) and (x, Az/2), G^{x, Az/2) is the tangential wall force per unit 
area in this layer, given by 

G-(a:,A^/2)= 1"^^ dzg^{x,z) 
Jo 

in terms of the tangential wall force density g^, and Gl{x, Az/2) is the tangential fluid force 
per unit area in the same layer, given by 

Gl{x, Az/2) = J^^ dz[d^a^^{x, z) + d^a^xix, z)] 

(6) 

^dx!^^ dza^^{x, z) + a^^{x, Az), 

where axx and azx are the xx and zx components of fluid stress in the dynamic state. Again, 
the fact that azx{x,0) — has been used in Eq. (6). According to Eq. (4), an integrated 
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form for the hydrodynamic tangential force balance may be obtained by subtracting Eq. (2) 
from Eq. (5): 

G:(a;,A;./2) + G{(a;,A;./2)=0, (7) 

In Fig. 9 we show the profiles of G"^ and —G^ for the boundary layer at the lower fluid-solid 
interface. It is seen that Eq. (7) is verified. 

Finally, we emphasize that the MD verification of Eqs. (2), (5), and (7) needs a rehable 
evaluation of the tangential fiuid forces G^, and G^.. It is essential to use the expressions 
in Eqs. (3) and (6), together with a stress measurement scheme that is reliable near the 
fiuid-solid interface (see Appendix C). 

C. Boundary Layer and Sharp Boundary Limit 

From Eqs. (3), (6), and (4), we have the hydrodynamic tangential fluid force per unit 
area 

Gl{x, /\z/2) = J^^ dz[d^a^^{x, z) + d^a^^{x, z)\ 

(8) 

= Jo dza^^{x, z) + a^^{x, Az), 

in the boundary layer at the lower fluid-solid interface. Integrated forms for Gl, and Gl 
are necessary because the tangential wall force is distributed within a finite distance from the 
wall. MD measurements show that, beyond the boundary layer (the first horizontal layer at 
the wall), the hydrodynamic tangential wall force G^ vanishes. That is, the tangential wall 
force density is distributed within the boundary layer only, and hence G^i^j Az/2) fully 
accounts for the hydrodynamic tangential wall force (at the lower fiuid-solid interface). The 
short-range nature of is due to the fact that as a hydrodynamic quantity defined by Eq. 
(4) , G^ solely arises from the "roughness" of wall potential through kinetic collisions of fiuid 
molecules with wall molecules [5]. The modulation amplitude for this roughness decreases 
quickly with increasing distance from the wall. Hence the hydrodynamic tangential wall 
force tends to saturate at a relatively short range (see Appendix B). Considering that G'^ 
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exists in the boundary layer only, hereafter we use G^i^) (— ^xi^^ ^^/'^)) denote the 
total hydrodynamic tangential wall force per unit area. 

We take the sharp boundary limit by letting strictly concentrate at ^ = 0: g^i^, z) = 
G'^{x)S{z) with the same G'^{x) per unit area. Rewriting Gl in Eq. (8) as 

Gax,Az/2) 
= J^^ dz[d^a^^{x, z) + d^a^^{x, z)] 
= ^zx{x, 0+) + /o+'^ dz[d^a^^{x, z) + d^d^^{x, z)\, 

we obtain 

G{(a:,A^/2) = a,,(a:,0+), 

because now local force balance requires d^^xx + dzdzx — above z — Q'^ . Hereafter we use 
Gi{,x) (= G^(a;, Az/2)) to denote the hydrodynamic tangential fluid force per unit area in 
the boundary layer. Therefore, 

Gi{x) = ~azx{xX) = -G':{x). (9) 

In the sharp boundary limit, the boundary layer extends from 2; = 0~ to 2; = O"*", and d^x 
varies from azx{x, 0^) = to cFzx{x, 0^) = G^(a;) at 2; = such that 

{V-a)-±^Gi5{zl 

in balance with the tangential wall force density g'^{x,z) = G'^{x)6{z). Equation. (9) may 
serve as a boundary condition in a continuum hydrodynamic model provided the continuum 
forms of d-zx{x, 0+) and are given. In Sec. IV A we show the validity of the Navier slip 

model for G^{x). In Sec. IV B we show the validity of the Newtonian constitutive relation 
for azx{x,0+). 

IV. TOWARDS A CONTINUUM HYDRODYNAMIC MODEL 

A. Local Slip Coefficient 

Hydrodynamic motion of fluid at the fluid-solid interface is represented by the slip ve- 
locity v^^*^. In MD simulations v^^*^ is obtained from the tangential fluid velocity in the 
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boundary layer, measured relative to the wall, i.e., v^'*^ = Vx at the lower fluid-solid inter- 
face or vf^P = Vx — V Sit the upper interface. The hydrodynamic viscous coupling between 
fluid and solid is described by a constitutive relation between the hydrodynamic tangential 
wall force (7^ and the slip velocity v^'''^^: 



where P{x) is the local shp coefficient at x. Equation (10) is generally referred to as the 
Navier slip model. The validity of this model has been verified by many MD studies, yet 
most of them have been done for homogeneous solid surfaces only [4-7]. In the present 
problem, the local nature of the Navier shp model can be better revealed. Obviously, the 
slip coefficient /3 must vary along the patterned surface where the fluid-solid interaction is 
modulated. In particular, being a local coupling constant, /3 should depend on x through 
the x-dependent parameter 5fs in potential Ufs, i.e.. 



where /3 as a function oi Sfs, denoted by P{Sfs), relates a microscopic interaction parameter 
and a coefficient in the hydrodynamic shp model. Equation (11) reflects the short-range 
nature of the fluid-solid interaction, and has been verified as follows. 

A series of independent MD simulations have been carried out, for Couette flow between 
two identical homogeneous solid surfaces, using two identical fluid-solid interaction potentials 
with a constant 6fs- Other simulation parameters remain unchanged, including p, T, p^, e/^, 
and afs- For each particular value of Sfg, the slip coefficient /? is obtained by measuring Gf^ 
and v^'''^^, each being a constant along the x direction. The functional dependence of /3 on 5/^, 
P{Sfs), has been numerically obtained by scanning a set of values for Sfg, as shown in Fig. 
10. (This dependence is qualitatively consistent with the large slip effect at a nonwetting 
fluid-solid interface [6].) Substituting into P{Sfs) the 5fs{x) profile in Eq. (1) then yields 
f3{x) in Eq. (11) for the patterned surface. Using this l3{x), plus the slip velocity profile 
yshp^^^ directly measured in the MD simulation (see Fig. 11), we calculate the hydrodynamic 



G:(x) = -(3(x)vf^(x), 



(10) 



P{x)^P{Sfs{x)), 



(11) 
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tangential wall force at the patterned surface according to Eq. (10). This calculated 
is shown in Fig. 12, and compared to the measured (which has already been given in Fig. 
7). The good agreement clearly verifies the local Navier slip model expressed by Eqs. (10) 
and (11). Of course, the local dependence of /3 on expressed by Eq. (11), is not without 
limit. Physically, (3 is defined as a phenomenological parameter over a microscopic length 
scale that is about a few lattice constants of the atomistic wall. Therefore, to validate Eq. 
(11) for I3{x) in Eq. (10), the lateral variation of Sfg must be slow. That is, the characteristic 
distance over which 5fs varies must be much larger compared to the wall lattice constant. 
That's why in Eq. (1) the modulation period P is chosen to be larger than lOcr. In fact, 
the disagreement between MD and continuum results has already been observed when the 
period of surface pattern approaches a molecular scale [17]. 



B. Newtonian Stress 

To obtain a continuum hydrodynamic model for fluid flow past a chemically patterned 
solid surface, we need a momentum transport equation and a slip boundary condition. Far 
away from the solid surface, the molecular density is uniform, the (hydrodynamic) stress 
is Newtonian, and the momentum transport can be well described by the Navier-Stokes 
equation 

^-Vp + r]V\, (12) 

where the viscosity rj = I.X^femj a"^ has been determined by MD measurements. Close to 
the rigid wall, however, the fluid density shows a short-range oscillation along the fluid-solid 
interface normal [z) across a few molecular layers [21]. It has been shown that, in the 
presence of this normal density variation, MD flow fields still can be well reproduced by 
continuum calculations using a constant viscosity, unless the channel is extremely narrow 
{E^^a) [22,23]. 

In the present problem, the fluid density shows a lateral oscillation (along x) as well, 
imposed by the patterned solid surface. With the lateral and normal density variations both 
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present near the solid surface, we proceed by using a constant viscosity throughout the fluid 
space in our continuum model. This is based on the observation that the magnitude of 
lateral density variation is on the same order of that of normal variation. To see the validity 
of the constant viscosity approximation, the hydrodynamic tangential stress azx{x, z) has 
been measured at 2; 2Az where density oscillation is f« 10%. This measured quantity is 
then compared to that calculated from the Newtonian relation rj{dzVx + dxVz) using the MD 
measured velocity held (see Fig. 13). The overall agreement supports the use of a constant 
viscosity. 

The hydrodynamic slip boundary condition is obtained by combining Eqs. (9), (10), 
and (11) with the Newtonian constitutive relation for the tangential stress a^x at the solid 
surface. This leads to the Navier shp boundary condition 

ms{x))vfP = -VdnVx, (13) 

with a local slip coefficient P{5fs{x)) (here n denotes the outward surface normal, —z at 
the lower surface or z at the upper surface). A local slip length, ls{x), can be defined by 
h{x)^rj/P{Sfs{x)). 

To summarize, our hydrodynamic model consists of the Navier-Stokes equation (12), 
the NBC (13), the impermeability boundary condition Vn = 0, and the incompressibility 
condition V • v = 0. Continuum calculations involve six parameters, including the system 
dimensions L along x and H along z, the shearing speed V, the ffuid density p, the viscosity 
7], and the slip coefficient f3 (as a function of Sfs)- We emphasize that our model has no 
adjustable parameter: the values for L, H, V, p come directly from MD simulations, 77 has 
been determined to be 2.1v^em/cr^, and (3 has been measured as a function of 5fs (see Fig. 
10). A continuous dependence of /3 on can be obtained from the MD data (in Fig. 10) 
through interpolation. 
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V. CONTINUUM HYDRODYNAMICS RESULTS 



A second order finite-difference scheme is designed to solve the Navier-Stokes equation. 
Essentially, it is a modified version of the pressure- Poisson formulation [24], where the in- 
compressibility condition is replaced by the pressure Poisson equation plus a divergence-free 
boundary condition for the velocity. The pressure-Poisson formulation has already been 
used to solve the Navier-Stokes equation for immiscible two-phase flow in the vicinity of 
the moving contact line, where a generalized NBC is involved [13,25]. (For details of the 
numerical algorithm, see Appendix C of Ref. [25]. It can be easily simplified to work for the 
present single-phase problem.) 

A. Comparison of MD and Continuum Results 

The two equations of our continuum model, Eqs. (12) and (13), come directly from the 
MD verification of the Newtonian stress (see Fig. 13) and the local NBC (see Eqs. (10) 
and (11), plus Fig. 12). Therefore, predictions from the continuum model are expected to 
match MD results quantitatively. 

Figure 14 shows the periodic profiles close to the lower wall of sinusoidal periodic 
fiuid-solid interaction. 

Figure 15 shows the periodic Vx profiles close to the lower wall of stepwise periodic fiuid- 
solid interaction. (This corresponds to the transverse case studied in Ref. [17].) In the first 
period -P/2 <x< P/2, Sfs{x) ^ 6o - 5i for -P/4 <x< P/4 or So + Si elsewhere. 

The overall agreement in the above two cases affirms the vahdity of the local NBC and 
the hydrodynamic model. The small discrepancies seen in Figs. 14 and 15 may be attributed 
to the relatively small H (= 13.6cr), the nonuniform boundary- layer fiuid, and the relatively 
fast variation of Sfg (especially in the sinusoidal case). 
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B. Adjustable Slip Length 



Consider a chemically patterned surface composed of a periodic array of stripes parallel 
to the y axis. Each stripe is of type A or S, and the periodicity along the x direction 
is realized by arranging the stripes according to • • -ABABAB ■ ■■ (see Fig. 1). This way 
of patterning has been modeled in MD simulations by a stepwise modulation of fluid-solid 
interaction, and the MD and continuum results agree well (see Fig. 15). 

The continuum model uses two different slip lengths Isa and Isb for the A and B stripes, 
respectively. While the chemical property of each stripe is kept invariant, the surface pattern 
may be continuously changed by varying the stripe widths. With the two slip lengths fixed, 
continuum flow fields have been calculated for stripe widths varied in a wide range. It is 
found that the slip behavior can continuously change from the regime of distinct slip domains 
(for large stripe widths) to the regime of effective medium (for small stripe widths, with an 
effective slip length showing up). Note that in the latter regime, the stripe widths, though 
small compared to the slip length, are still considered large enough to justify the vahdity of 
the continuum model. The effective-medium regime provides a way to systematically tune 
the effective slip length. A desired slip length, if not naturally produced at a chemically 
homogeneous surface, can thus be realized at a patterned surface. 

Continuum hydrodynamic calculations have been carried out for Couette flow between 
two solid planar surfaces parallel to the xy plane. The lower surface is patterned and 
the upper surface is homogeneous. The fluid is sheared by flxing the lower surface and 
translating the upper surface at a constant velocity V in the x direction. We focus on the 
regime of viscous flow with small Reynolds numbers {Re = pVH/r] pa 1, where H is the 
distance between the two solid surfaces) . With the fluid velocity v measured by V and the 
coordinates x and z measured by H, the scaled steady-state velocity flelds, v{x/H, z/H)/V, 
are controlled by the dimensionless parameters Isa/H, Isb/H, wa/H, and wb/H, where wa 
and wb are the widths of A and B stripes, respectively. 

Figures 16 and 17 show the steady-state shp proflles at six differently patterned surfaces. 
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obtained for Isa/H — 0.147 and Isb/H — 0.441. The values of wa/H and wb/H have 
been varied to show the control of an effective slip length by surface patterning. In the 
present Couette-flow geometry, the slip velocity has to oscillate in the x direction at the 
lower surface, where the slip length is modeled as a periodic function of x according to the 
surface patterning. Prom Fig. 16, it is obvious that the shp velocity becomes less oscillatory 
in magnitude when the pattern period P is decreased. In fact, for P/i7 = 1/4 and 1/8, the 
relative magnitude of oscillation for Vx, defined as {v'^^^ — v'^^'")/2v^^^, becomes 3.3% and 
1.5%, respectively, while for P/H — 1, that oscillation magnitude is 14%. Therefore, when 
the solid surface is patterned with a sufficiently small period, it appears to be homogeneous. 
The nearly uniform shp velocity can be used to define an effective shp length l^^f. For the 
present Couette fiow, it is given by 

H ~ V - v"^"" ' 

which yields 1%^^ jB. — 0.196 from the data in Fig. 16. For sufficiently small pattern period 
P, the effective slip length can be further tuned by varying the ratio of the A stripes (with 
a small slip length) to the B stripes (with a large shp length). Figure 17 shows that by 
fixing P/ H Ski and increasing the proportion of the A stripes, the average amount of slip 
is appreciably reduced. That is, the effective slip length l^^^ decreases with the increasing 
proportion of the A stripes. 

VI. CONCLUSION AND DISCUSSION 

In this paper we have used both the MD and continuum simulations to study the slip 
boundary condition for single-phase fiow past a chemically patterned surface. Modulation 
of fiuid-solid interaction along the patterned surface induces a lateral variation in various 
stresses and forces in the fiuid. This has been observed in both the static equilibrium state 
and the dynamic steady state. Given this lateral variation, a subtraction scheme has been 
developed for extracting hydrodynamic quantities from MD simulations. Based on this 
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scheme, the vahdity of the Navier boundary condition has been verified in describing the 
local hydrodynamic viscous coupling between the fiuid and the patterned solid surface. A 
continuum hydrodynamic model has been formulated. It is capable of yielding steady-state 
fiow fields in quantitative agreement with those from MD simulations. We have used the 
continuum model to show that an effective slip length can be realized and tuned by surface 
patterning. 

The framework of heterogeneous multiscale method [26] provides a general scheme for 
developing hybrid atomistic-continuum methods. The hybrid methods use MD to extract 
boundary condition(s) and/or constitutive relation(s) for continuum calculation. Whatever 
a continuum calculation needs, the relevant hydrodynamic quantities have to be properly 
extracted from MD data. In the present problem, we see that a subtraction of the static 
component from the dynamic quantity is essential to obtaining the correct hydrodynamic 
information. We emphasize that hybrid methods are mostly used to investigate problems 
involving complex fiuids, complex fiuid-solid interactions, and/or complex geometries, while 
in such problems, the static configurations are usually nontrivial in a molecular point of 
view. 
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APPENDIX A: ATOMISTIC SOLID WALLS 

The MD results presented in this paper were obtained from simulations using solid walls 
constructed by two [001] planes of an fee lattice. We have also carried out MD simulations 
using thicker confining walls. The number of crystalline layers ([001] planes of fee lattice) is 
changed from two to four in constructing each of the two walls. The fiuid-solid interaction 
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potential Ufs is then cut off at a larger distance Tc = 2.5a. It turned out that both the 
velocity field and the stress field remain essentially the same. 

APPENDIX B: TANGENTIAL WALL FORCE 

In terms of the tangential wall force density g'^, we define the tangential wall force per 
unit area in a horizontal layer as G'^{x,^{zi + Z2)) = J^'^ dzg'^{x, z), which is the total 
tangential wall force accumulated across that layer, sandwiched between z = zi and z = Z2. 

The wall force can be singled out by separating the force on each fluid molecule into 
fluid-solid and fluid-fluid components. The fluid molecules close to the crystalline wall 
experience a fast periodic modulation in interaction with the wall. When coupled with 
kinetic collisions with the wall molecules, there arises the first component in the tangential 
wall force, sharply distributed along z close to the wall. In the presence of a periodically 
varied fluid-solid interaction (at the lower fluid-solid interface, with a period much larger 
than that of the crystalline structure), the second component in the tangential wall force 
arises along with a slow oscillation of fluid molecules' distribution. This component is also 
narrowly distributed along z near the wall, but with a distribution width larger than that 
of the flrst component. This may be understood as follows. 

Two quite different periods exist in the periodic modulation in fluid molecules' inter- 
action with the wall. The first period is from the crystalline structure of the wall. For a 
fiuid molecule close to the sohd wall, the interaction with one particular (the closest) wall 
molecule can be much stronger than all the others. As this fiuid molecule moves laterally 
but remaining close to the wall, it would thus experience a strong periodic modulation in its 
interaction with the wall. This lateral inhomogeneity is generally referred to as the "rough- 
ness" of the wall potential [5] . Away from the fiuid-sohd interface, each fiuid molecule can 
interact with many wall molecules on a nearly equal basis. Thus the modulation amplitude 
of the wall potential would clearly decrease with increasing distance from the wall. Hence 
the first component of the tangential wall force tends to saturate at a relatively short range. 
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The second period, being much larger than the first one, comes directly from the fiuid-solid 
interaction potential. As a consequence, the distribution of fluid molecules would be pe- 
riodically adjusted along the wall. This occurs within a relatively long distance from the 
surface, as determined by the range of fluid-solid interaction and the correlation length of 
fluid. Therefore, the second component of the tangential wall force has a wider distribution 
than the first one. 



APPENDIX C: STRESS MEASUREMENT 

The Irving-Kirkwood expression [27] has been widely used for stress measurement in 
MD simulations. However, as pointed out by the authors themselves, the leading-order 
expression for the interaction contribution to the stress tensor represents only the leading 
term in an asymptotic expansion, accurate when the interaction range is small compared 
to the range of hydrodynamic variation. As a consequence, this leading-order expression is 
not reliable near a fiuid-fiuid or a fiuid-solid interface. In this paper, the stress components 
Uxx and Uzx are obtained from the time averages of the kinetic momentum transfer plus 
the fiuid-fiuid interaction forces across the fixed x and z bin surfaces. In particular, we 
directly measure the x component of the fiuid-fiuid interaction forces acting across the x{z) 
bin surfaces, in order to obtain the xx{zx^ component of cr^j, the interaction contribution to 
the stress tensor a. Here "acting across" means that the line connecting a pair of molecules 
intersects the bin surface (the so-called Irving-Kirkwood convention [27]). For more details, 
see Appendix B in Ref. [25]. 
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FIGURES 



o 




FIG. 1. Schematic of the sohd planar surface which is parallel to the xy plane and patterned 
along the x direction. Interaction between a fluid molecule (empty circle) and a solid molecule 
(solid circle) separated by a distance r is modeled by the potential Ufs{r), with the parameter 6fs 
varying with the equilibrium x coordinate of the solid molecule. 
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FIG. 2. Schematic of the Couette-flow geometry adopted in MD simulations. The lower solid 
surface is patterned in the x direction and the upper surface is homogeneous. Moving the upper 
wall with speed V along x generates the Couette flow. 
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FIG. 3. Static tangential wall force plotted as a function of x. The circles de- 

note G^^ = dzg^^{x,z) measured in the boundary layer at the lower fluid-solid interface, 
showing a fast oscillation superimposed on a slow variation of period P; the squares denote 
Gx'^ = JH-Azd-Z9x^{x, z) measured in the boundary layer at the upper interface, showing a fast 
oscillation only. Here g'^^ is the tangential wall force density. 
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FIG. 4. z-intcgrated static normal stress across the boundary layer, S^^ = dza^^{x, z), 
plotted as a function of x. 
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FIG. 5. Static tangential stress a'^j.{x,z) at z = Az, plotted as a function of x. 
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FIG. 6. Static tangential wall force G'^^{x, Az/2) and the negative of static tangential fluid 
force —GIP{x, Az/2) measured in the boundary layer at the lower fluid-solid interface. The circles 
denote (see Fig. 3); the squares denote It is seen that + G^P = (Eq. (2)). Our 

stress measurement scheme has smoothed the kinetic contribution to a^^. That's why the fast 
oscillation in — is less prominent than that in G^^. 
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FIG. 7. Hydro dynamic tangential wall force G^(x,Az/2) (circles) and stress azx{x, Az) 
(squares) for the boundary layer at the lower fluid-solid interface. Note that the magnitude for the 
variation in each hydrodynamic quantity is much smaller than that seen from the corresponding 
static quantity (see Figs. 3 and 5). 
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FIG. 8. z-integrated hydro dynamic normal stress across the boundary layer at the lower 
fluid-solid interface, S^^. = dza^^{x, z), plotted as a function of x. Note that the magni- 
tude for the variation in this hydrodynamic quantity is much smaller than that seen from the 
corresponding static quantity (see Fig. 4). 
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FIG. 9. Hydrodynamic tangential wall force G^(a;, Az/2) and the negative of hydrodynamic 
tangential fluid force — G|(x,A2;/2) in the boundary layer at the lower fluid-solid interface. The 
circles denote G'^ (see Fig. 7); the squares denote — G^. It is seen that G^ + G^ = (Eq. (7)). 
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FIG. 10. Slip coefficient (3 plotted as a function of the parameter Sfg in Ufs{r). 
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FIG. 11. Slip velocity v^^'^^{x) at the patterned surface, plotted as a function of x. It is seen 
that v^^P varies periodically along x, being small in the large 5fs (large (3) region and large in the 
small 5fs (small (3) region. Larger the (3 coefHcient, smaller the slip length and slip velocity. The 
no-slip boundary condition corresponds to the limit of /? — >^ oo. 
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FIG. 12. Hydrodynamic tangential wall force G'^(x) at the lower fluid-solid interface. The 
circles denote G^{x) directly measured in the MD simulation (see Fig. 7); the squares denote 
G'^{x) calculated from P{x) and the measured vf'^'P{x) according to Eq. (10). 
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FIG. 13. Hydrodynamic tangential stress azx{x, z) at z = 2Az, plotted as a function of x. The 
circles denote azx{x, 2Az)—azx{x, 2Az), obtained according to the subtraction scheme for MD data; 
the dashed line denotes r]{dzVx + dxVz), calculated from the MD flow field. Here rj = 2.3-v/em/(T^ 
has been used to achieve the best agreement. 
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FIG. 14. Vx plotted as a function of x for four z levels close to the lower wall of sinusoidal 
fluid-solid interaction (see Eq. (1) for 5fs{x) of period P). The symbols denote the MD data and 
the lines represent the corresponding continuum results. Here the MD data are obtained for a 
system of H = 13.6(7, L = 2P = 54.4(7, and V = 0.25 y/ejm. The four z levels are at ^ = 0.425(7 
(circles and solid line), 1.275(7 (squares and dashed line), 2.125(7 (diamonds and dotted line), and 
2.975(7 (triangles and dash-dotted line). 
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FIG. 15. Vx plotted as a function of x for four z levels close to the lower wall of stepwise 
fluid-solid interaction. The symbols denote the MD data and the lines represent the corresponding 
continuum results. Here the MD data are obtained for a system of H = 13.6cr, L = P = MAa, 
and V = 0.25^e/m. The four z levels are at z = 0.425(7 (circles and solid line), 1.275(7 (squares 
and dashed line), 2.125(7 (diamonds and dotted line), and 2.975(7 (triangles and dash-dotted line). 
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FIG. 16. Slip profiles at differently patterned solid surfaces. Here the scaled slip velocity Vx/V 
is plotted as a function of x/H. The pattern period P = wa + wb varies from H to H/8, with 
equal stripe widths {wa = wb)- The solid line is ioi P/H = 1, the dashed line for P/i? = 1/2, the 
dotted line for P/H = 1/4, and the dash-dotted line for P/H = 1/8. 
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FIG. 17. Slip profiles at differently patterned solid surfaces. Here the scaled slip velocity Vx/V 
is plotted as a function of x/H. The pattern period P = wa+wb equals H/4, with wa/wb varying 
from 1/3 to 3. The solid line is for wa/wb = 1/3, the dashed line for wa/wb = 1 (same as the 
dotted line in Fig. 16), and the dotted line for wa/wb = 3. 
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